package provider.impl.coalescent.me;

import coalescent.data.K69Data;
import coalescent.model.K69;
import commons.core.ModelEstimator;
import commons.core.Singletons;

import java.math.BigDecimal;

import static commons.core.Singletons.MATHCONTEXT_128_HALF_UP;

/**
 * @author Susanta Tewari
 * @version version
 * @history Created on 2/10/13
 * @since 1.4.0
 */
final class K69Estimator_Watterson
        extends ModelEstimator<K69, K69Data, K69Data, K69Estimator_Watterson> {

    public K69Estimator_Watterson(final K69Data statistic) {
        super(statistic);
    }

    @Override
    public K69 getEstimate() {

        final K69Data data      = getData();
        final int mutationCount = data.getMutationCount();
        final int n             = data.getSampleSize();
        final double estimate = new BigDecimal(mutationCount).divide(harmonic_n(n - 1),
                                    MATHCONTEXT_128_HALF_UP).doubleValue();

        return new K69(estimate);
    }

    /**
     * @return true
     */
    @Override
    public boolean isVarianceKnown() {
        return true;
    }

    @Override
    public BigDecimal getVariance() {

        final BigDecimal theta = new BigDecimal(getEstimate().getMutationRate());
        final int n            = getData().getSampleSize();
        final BigDecimal h1    = harmonic_n(n - 1);
        final BigDecimal h2    = harmonic_n2(n - 1);
        final BigDecimal f1    = theta.divide(h1, MATHCONTEXT_128_HALF_UP);
        final BigDecimal f2 = BigDecimal.ONE.add(theta.multiply(h2).divide(h1,
                                  MATHCONTEXT_128_HALF_UP));

        return f1.multiply(f2);
    }

    static BigDecimal harmonic_n(final int n) {

        BigDecimal result = BigDecimal.ONE;

        for (int i = 2; i <= n; i++) {

            result = result.add(BigDecimal.ONE.divide(new BigDecimal(i),
                    Singletons.MATHCONTEXT_128_HALF_UP));
        }

        return result;
    }

    static BigDecimal harmonic_n2(final int n) {

        BigDecimal result = BigDecimal.ONE;

        for (int i = 2; i <= n; i++) {

            result = result.add(BigDecimal.ONE.divide(new BigDecimal(i * i),
                    Singletons.MATHCONTEXT_128_HALF_UP));
        }

        return result;
    }

    @Override
    public K69Estimator_Watterson factory(final K69Data data) {
        return new K69Estimator_Watterson(data);
    }
}
